home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libfft / TRY / c_check1d.c next >
Encoding:
C/C++ Source or Header  |  1994-08-02  |  5.1 KB  |  236 lines
  1. #include <stdio.h>
  2. #include <sys/time.h>
  3.  
  4. #include "fft.h"
  5. #include "constant.h"
  6.  
  7. /*                        */
  8. /*    Precision dependant declarations    */
  9. /*                        */
  10. #ifdef    ZOMPLEX
  11.  
  12. #define TOLERANCE   DTOLERANCE
  13. typedef        double    this_half;
  14. typedef        zomplex    this_type;
  15.  
  16. #define        THIS_SUM    dsum_
  17. #define        THIS_GEN    zgen_
  18. #define        THIS_FT    zft1d_
  19.  
  20. #define        THIS_FFTI    zfft1di
  21. #define        THIS_FFT    zfft1d
  22. #define        THIS_SCAL    zscal1d
  23.  
  24. #endif
  25. #ifdef    COMPLEX
  26.  
  27. typedef        float        this_half;
  28. typedef        complex        this_type;
  29.  
  30. #define TOLERANCE   STOLERANCE
  31.  
  32. #define        THIS_SUM    ssum_
  33. #define        THIS_GEN    cgen_
  34. #define        THIS_FT    cft1d_
  35.  
  36. #define        THIS_FFTI    cfft1di
  37. #define        THIS_FFT    cfft1d
  38. #define        THIS_SCAL    cscal1d
  39.  
  40. #endif
  41.  
  42. /*                        */
  43. this_half THIS_SUM();
  44.  
  45. void inimat_();
  46. void GetArguments();
  47. void get_values();
  48.  
  49. int size, inc, n_trials;
  50. this_type *pa, *pb, *pRef, *pWork;
  51.  
  52. main(argc,argv)
  53. int argc;
  54. char *argv[];
  55. {
  56.     int i, j, k, n_total, is_wrong, nerr;
  57.     double x, y, dx, dy, emax, sum, err, maxerr;
  58.     this_half    t;
  59.  
  60. /* ******************************************************* */
  61.     GetArguments( argc, argv);
  62. /* ******************************************************* */
  63.  
  64.     srandom( (123*getpid()) | 0x01);
  65.  
  66.     for( k = 0; k < n_trials ; k++) { 
  67.     get_values();
  68.  
  69.     if( !(k%1) && (k)) printf("\n");
  70.     printf( "<%3d,%3d>", size, inc);
  71.     fflush(stdout);
  72.  
  73.     n_total = ((size+1)*inc);
  74.     pa = (this_type *)malloc( (3 * n_total + size + FACTOR_SPACE) * sizeof(this_type));
  75.     if( !(pa) ) {
  76.         fprintf( stderr, "Could not allocate ... Exiting\n");
  77.         exit (-1);
  78.     }
  79.     pb = pa + n_total;
  80.     pRef = pb + n_total;
  81.     pWork = pRef + n_total;
  82.  
  83. /* *******************************************************
  84.     Initialize arrays
  85. ******************************************************* */
  86.     THIS_GEN(pRef, &n_total);
  87.     bcopy( pRef, pa, n_total*sizeof(this_type));
  88.     bcopy( pRef, pb, n_total*sizeof(this_type));
  89.  
  90. /* *******************************************************
  91.     direct Fourier Transform
  92. ******************************************************* */
  93.         j = -1;
  94.     THIS_FT( &j, &size, pb, &inc);
  95.     pWork = THIS_FFTI( size, pWork);
  96.     is_wrong = THIS_FFT( -1, size, pa, inc, pWork);
  97.  
  98.     emax = TOLERANCE*sqrt(size);
  99.     for(i = 0, nerr=0 ; i < (size*inc) ; i ++) {
  100.         x = pa[i].re - pb[i].re;
  101.         y = pa[i].im - pb[i].im;
  102.         if( (t= x*x + y*y) > (emax)) {
  103.         fprintf( stderr, "@", size, i, t);
  104.         nerr++;
  105.         break;
  106.         }
  107.     }
  108.     if( !(nerr)){
  109.         fprintf( stderr, ".", size, i, t);
  110.     }
  111.     i = inc * 2;
  112.     x =  THIS_SUM( &size, pRef, &i);
  113.     y =  THIS_SUM( &size, ((char *)pRef) + sizeof(this_half), &i);
  114.     dx = x - pa[0].re;
  115.     dy = y - pa[0].im;
  116.     if( (t= dx *dx + dy*dy) > TOLERANCE){
  117.         fprintf( stderr, 
  118.         "\nError direct FFT(%d) at index #0 : (%f,%f)<->(%f,%f)error(%f,%f)",
  119.         size, pa[0].re, pa[0].im,x,y, dx, dy);
  120.  
  121.     }
  122.     dx = x - pb[0].re;
  123.     dy = y - pb[0].im;
  124.     if( (t= dx *dx + dy*dy) > TOLERANCE){
  125.         fprintf( stderr, 
  126.         "\nError direct DFT(%d) at index #0 : (%f,%f)<->(%f,%f)error(%f,%f)\n",
  127.         size, pa[0].re, pa[0].im,x,y, dx, dy);
  128.  
  129.     }
  130. /* *******************************************************
  131.     inverse Fourier Transform
  132. ******************************************************* */
  133.     is_wrong = THIS_FFT( 1, size, pa, inc, pWork);
  134.  
  135.     i = inc * 2;
  136.     x =  THIS_SUM( &size, pb, &i);
  137.     y =  THIS_SUM( &size, ((char *)pb) + sizeof(this_half), &i);
  138.     dx = x - pa[0].re;
  139.     dy = y - pa[0].im;
  140.     nerr = 0;
  141.     emax = size * TOLERANCE;
  142.     if( (t= dx *dx + dy*dy) > TOLERANCE){
  143.         fprintf( stderr, "#");
  144.         nerr ++;
  145.     }
  146.     if( !(nerr)){
  147.         fprintf( stderr, ".", size, i, t);
  148.     }
  149.  
  150.     t = 1. / (double)size;
  151.     THIS_SCAL( size, t, pa, inc);
  152.  
  153.     emax = TOLERANCE;
  154.     emax = emax * emax;
  155.  
  156.     sum = 0.;
  157.     err= 0.;
  158.     nerr= 0;
  159.     maxerr= 0.;
  160.     for(i = 0 ; i < (size*inc) ; i += inc) {
  161.         sum += (pRef[i].re * pRef[i].re) + (pRef[i].im * pRef[i].im);
  162.         x = pa[i].re - pRef[i].re;
  163.         y = pa[i].im - pRef[i].im;
  164.         if( (t= x*x + y*y) > (emax))
  165.         nerr++;
  166.         err += t;
  167.         if( t > maxerr) maxerr = t;
  168.     }
  169.     err = sqrt(err / (double)size);
  170.     sum = sqrt(sum / (double)size);
  171.     maxerr = sqrt(maxerr);
  172.     printf(": avg:(%8.3g /%8.3g)= %8.3g <> maxerr: (%8.3g /%8.3g)= %8.3g", 
  173.     err, sum, err/sum, maxerr, sum, maxerr/sum);
  174.     if(nerr){
  175.         printf("\n%d errors detected \n");
  176.         exit(-2);
  177.     }
  178.     free(pa);
  179.     }
  180.     printf("\n");
  181.     return(0);
  182. }
  183.  
  184. int is_random;
  185.  
  186. void GetArguments( argc, argv)
  187. int argc;
  188. char *argv[];
  189. {
  190.     int i, j, k;
  191.     int nerror = 0;
  192.  
  193. #define ON    1
  194.  
  195.     n_trials = DEF_TRIALS;
  196.     is_random = 1;
  197.  
  198. /* ******************************************************* */
  199.     for ( i = 1 ; (i < argc) && (nerror != ON) ; i ++ ) {
  200.     if( argv[i][0] == '-') {
  201.         switch ( argv[i][1]) {
  202.         case 'i' :
  203.         case 'I' :  
  204.                 is_random = 0;
  205.                 break;
  206.         default  :  nerror = ON;
  207.         }
  208.     }
  209.     else {
  210.         if( i+1 > argc)
  211.         nerror = ON;
  212.         else { 
  213.         n_trials = atoi( argv[i]);
  214.         }
  215.     }
  216.     }
  217.     if( nerror == ON) {
  218.     fprintf( stderr, 
  219.         "Usage : %s [-p(arallel)] [n_trials]\n", argv[0]);
  220.     exit(-1);
  221.     }
  222. }
  223.  
  224. void get_values()
  225. {
  226.   if( is_random){ 
  227.     size = random() % MAX_SIZE + 1;
  228.     inc = random() % MAX_STRIDE + 1;
  229.   } else{
  230.     printf( "Enter SIZE : ");
  231.     scanf("%d", &size);
  232.     printf( "Enter STRIDE : ");
  233.     scanf("%d", &inc);
  234.   }
  235. }
  236.